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SUMMARY 

A second-order finite difference and two spectral methods including a Cheby- 
shev tau and a Chebyshev collocation method have been implemented to determine 
the linear hydrodynamic stability of an unbounded shear flow. The velocity profile 
of the basic flow in the stability analysis mimicks that of a one— stream free mixing 
layer. Local and global eigenvalue solution methods are used to determine individ- 
ual eigenvalues and the eigenvalue spectrum, respectively. The calculated eigen- 
value spectrum includes a discrete mode, a continuous spectrum associated with 
the equation singularity and a continuous spectrum associated with the domain 
unboundedness. The efficiency and the accuracy of these discretization methods in 
the prediction of the eigensolutions of the discrete mode have been evaluated by 
comparison with a conventional shooting procedure. Their capabilities in mapping 
out the continuous eigenvalue spectra are also discussed. 
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1. INTRODUCTION 


This paper studies the numerical solutions of a boundary value problem using 
four different solution methods. The Rayleigh equation 1 governs the inviscid in- 
stability properties of linearized disturbances. With the conventional normal mode 
representations for disturbances, the equation for the complex amplitude, v, of the 
velocity perturbation in the y-direction is 

{( w - a 2 ) _ J o = 0 (1) 

In a spatial analysis, it is assumed that the disturbances with real frequency, u>, 
traveling at the speed, u?/cti r , are amplified at the rate exp(— a{x) upon the basic 
parallel flow described by the mean velocity in the ^-direction , U(y). a denotes the 
complex wavenumber and oj the frequency. The Rayleigh equation, together with 
the boundary conditions 


v — > 0 , y — > ± oo (2) 

defines the basic linear inviscid instability problem, in the form of a boundary value 
problem for parallel free shear flows in an unbounded domain. Thus we are solving 
an eigenvalue problem to determine the dispersion relationship 

a = a ( tx> ) (3) 

Traditionally, the hydrodynamic stability problem has been solved by shooting 
techniques 2 . This involves the solution of two initial value problems with the two 
boundary conditions as their, respective initial conditions. Eigenvalues are deter- 
mined by satisfying a matching condition at a certain intermediate point or bound- 
ary in the flow. A good knowledge of the characteristics of the solutions is often 
needed in the shooting technique in order to make a good initial guess for the eigen- 
values. Also, the accuracy of the solutions is often limited by the accuracy of the 
numerical integration scheme. 

Recently, interest in the application of spectral approximation methods has 
grown in all branches of science and engineering. Spectral methods are useful 
in problems where high resolution is required 3,4 . These methods have also been 
studied in classical problems of computational fluid dynamics 5,6 and in turbulence 
simulations 7 . Canuto et al. 8 have given a detailed description of the technical as- 
pects of the applications of various spectral methods in fluid dynamics. 

Another application of spectral methods has been in the solution of hydrody- 
namic stability problems. For example, Orszag 9 compared spectral methods based 
on different expansion polynomials and used a Chebyshev polynomial expansion to 
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study the temporal instability of plane Poiseuille flow. When this discretization is 
applied to the linearized equations of hydrodynamic stability with the appropriate 
boundary conditions, an algebraic eigenvalue problem is obtained. For temporal 
instability, in which a fixed wavelength disturbance grows or decays in time, the 
complex frequency is the eigenvalue. This parameter appears linearly in the prob- 
lem and standard algebraic eigenvalue techniques may be used to determine the 
eigenvalues. However, most shear driven instabilities are spatially unstable. In this 
case a disturbance of fixed real frequency grows or decays in space, and the com- 
plex wavenumber of the disturbance is the eigenvalue. The wavenumber appears 
nonlinearly in the problem and standard algebraic eigenvalue techniques are not ap- 
plicable. Bridges and Morris 10 showed that the Linear Companion Matrix Method 
(LCM) and a method based on matrix factorization (MF) could be applied success- 
fully to this problem. The LCM approximates the entire eigenvalue spectrum. The 
factorization method gives only a subset of the eigenvalue spectrum. However, the 
size of the companion matrix in the LCM is p times that of the original matrices 
where p is the order to which the eigenvalue appears. Note that neither of the two 
methods requires initial guesses for eigenvalue calculations. 

So far the applications of the spectral approximations in conjunction with the 
global eigenvalue solution methods have been limited to the boundary layer or 
bounded shear flow instability problem 10,11 . In the present analysis, two spectral 
methods, the Chebyshev tau method and a Chebyshev collocation^ method are used 
to determine the spatial, inviscid stability .of a free shear layer. The resulting 
eigenvalue problem in which the parameter appears nonlinearly is solved using the 
LCM and MF. 

Moreover, the advent of these computational methods, LCM and MF, has 
made it feasible to solve the spatial stability problem using other discretization 
techniques, such as finite difference formulation. The application of the finite dif- 
ference techniques also results in an eigenvalue problem in which the parameter 
appears nonlinearly. Therefore, standard algebraic eigenvalue techniques are not 
applicable in this case as well and the global eigenvalue solution methods have to 
be used as in the case of spectral approximations. For both spectral and finite dif- 
ference discretizations, the accuracy of the solutions depends mainly on the order 
of approximations. In the present analysis, the order of approximation is measured 
by either the number of Chebyshev polynomials used in the spectral methods or the 
number of grid points at which differential equations are discretized in the finite dif- 
ference formulation. Finite difference discretizations, however, are much simpler to 
implement than spectral methods. Therefore, the finite difference approximations 
provide a viable alternative for hydrodynamic stability analyses using the global 
eigenvalue solution methods. 

In the present paper two spectral approximations are used to determine the 
spatial, inviscid stability of a two-dimensional free shear layer. A simple finite 
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difference scheme is also considered. The solutions are compared with that us- 
ing the shooting procedure. This stability analysis is of interest as experimental 
observations 12 ’ 13 have shown that the local characteristics of large-scale coherent, 
turbulent structures in free mixing layers are described remarkably well as inviscid 
instability waves. If these observations are to form the basis of a turbulence model 
it is valuable to have efficient numerical schemes to solve the inviscid hydrodynamic 
stability problem. Liou 14 has successfully implemented these global approximation 
schemes in developing turbulence models based on a linear theory to simulate the 
evolution of a turbulent free mixing layer. 

In the following sections, the basic boundary value problem, the numerical 
discretizations and the eigenvalue solutions are first described. Comparisons of the 
accuracy and efficiency of the schemes are then given. The various features of 
the eigenvalue spectrum of the Rayleigh equation unveiled by the global n um erical 
approximations are also discussed. 

2. FORMULATION 

A transformation 

* = f (V ) (4) 

which maps the unbounded physical domain onto the finite Chebyshev domain [-1,1] 
must be used to apply the Chebyshev spectral methods. The transformed Rayleigh 
equation becomes 

[ Uv]a 3 — [u;f)]a 2 — [ U ( v'm )' m — ( U'm )' mv ] a+ iv ( v'rn )' m = 0 (5) 

where m denotes the metric .of the transformation and ( )' denotes d/d z. The 
boundary conditions become 


v -» 0 , z -*■ ± 1 (6) 

For the Rayleigh equation, it is shown below that the system of equations generated 
by each of the three approximation methods forms an eigenvalue problem with the 
eigenvalue, cc, appearing nonlinearly. That is, 

D 3 (a)v = 0 (7) 

where 

(a;) = Cock 3 + Cia 2 + C20; -j- C3. 

The C 0 ,C 1 ,C 2 and C 3 are the coefficient matrices of the lambda matrix D 3 (a). 
The components of the eigenvector v are either the expansion coefficients of the 
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Ch.ebysh.ev series approximation or the solution vectors themselves. The eigenvalues 
of the system are the roots of the characteristic equation 

det D 3 (a) = 0 (8) 


2.1 Chebyshev Tau Method 

In the Chebyshev tau (CT) method the solution v is approximated by a trun- 
cated finite series expansion of Chebyshev polynomials, 

v(z) ~ y(z) = ^ T »( 2 ) ( 9 ) 

n= 1 

where T n (^) is the n th order Chebyshev polynomial of the first kind. The vari- 
ous properties of Chebyshev polynomials can be found in Fox and Parker 15 . For 
convenience the Rayleigh equation is now written in integral form: 

Aa 3 + Boi 2 -f Coi + D -f- byz 4- 6 2 = 0 (10) 


where 



A = 

J J U v dzdz , 


B = 

— xv J J v dzdz, 

C — ~m 2 Uv 

3 

+ -z) 

^ \m 2 )'Uv dz + 2 J m 2 U'v dz— 

//£ 

)"Uv dzdz — j J (m 2 )'U'v dzdz, 

D — w {m 2 v 

-1/ 

{m 2 )'v dz + J f(^)%}dzdz, 


and b\ , b^ are integration constants. The series representation of v(z) is substi- 
tuted into equation (10) and the integrations are performed by making use of the 
Chebyshev relations 



r.M 

^n+l(z) 

2(n+l) 


T2{Z ) ] 

Tn-l(z) 


n — 0 
n = 1 

n > 2 


(ii) 
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After the application of the above expressions to the integrated Rayleigh equation, 
equation (10), we obtain (iV — 1) equations by equating the coefficients of the 
polynomials of degree n, n = 2, ..N. The resulting system of equations can 

be expressed in the form of equation (7). The details of the application of the 
Chebyshev tau method to hydrodynamic stability problems can be found in Bridges 
and Morris 10 and Liou 16 . 

2.2 Chebyshev Collocation Method 


It is well known that a smooth function, F(z), can be approximated by poly- 
nomials in 2 . The resulting polynomials, such as the Lagrange interpolation poly- 
nomial based on equally spaced collocation points, however, typically diverge as the 
number of collocation points increases. The poor convergence behavior of polyno- 
mial interpolation can be avoided by relating the collocation points to the structure 
of orthogonal polynomials, like Chebyshev or Legendre polynomials. In the Cheby- 
shev collocation (CC) method the collocation points, zj are the extrema of the iVth 
order Chebyshev polynomials Tjv(z). That is, 

Zj = cos(f), j = 0,...,iV (12) 

There are other choices of the collocation points 17 . The approximated solution 

becomes 

N 

v (z) = J2 v ( z i) b j( z ) (13) 

i=o 

where bj(z) axe the expansion orthogonal functions. The approximation simultane- 
ously interpolates the solution at each collocation point. That is, 

bj( z i) = $ij (14) 


The details of these expressions can be found in Voigt 17 . The resulting system 
of equations obtained by evaluating the differential equation, equation (5), at the 
collocation points can be put into the form of equation (7). 


2.3 Finite Difference Method 

The Rayleigh equation can also be discretized by finite difference (FD) methods 
The discretization is performed in the transformed plane where z G [-1,1]- The finite 
difference approximation to the Rayleigh equation at grid point i is 


~Ui di]a 3 + [uvi]a 2 + [/i( 


^8+1 Jpi— 1 '\ I X /^*+l ^^8 T ^8—1 \ I X * 1 

■ ) + M — ) + hVi]a 


2AS 


Az 2 
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(15) 


where 


. r r /^i+1 Vi— 1 \ . r / ^i+1 2 Vi -f- Vi— i 

+ i M 7T7TZ ) + M — 


2Az 


Az 2 


)} = 0 


fi = {U m m')i , f 2 = ( Um 2 )i , / 3 

/ 4 = (-wm 2 )i , / 5 = (-m 


= (-um m')j, 
(* = 1...N) 


and N denotes the total number of the grid points. The first and the second deriva- 
tives are approximated by corresponding second order finite difference formulas. 
It is important to note that there are no additional computational difficulties if 
higher order difference formulas are used. The finite difference discretization has 
been found to be more straight forward to formulate than either the Ghebyshev 
tau or the Chebyshev collocation methods described above. The application of the 
equation (15) at each grid point gives rise to a system of equations in the form of 
equation (7). 

In the Chebyshev tau method, the eigenvectors of the eigenvalue problem give 
the spectrum of the expansion. While in the collocation and the finite difference 
method, the eigenvectors are the solution vectors themselves. 


2.4 Boundary Conditions 


The boundary conditions for the present problem are 

v(±oo) = 0 (16) 

since the spatial instability waves must decay far from the shear layer. In the 
Chebyshev tau method, the boundary conditions become 

N 

y(±l) = f + Y, (i 1 )"”" = 0 ( 17 ) 

n— 1 

The addition of these two equations to the set of equations formed in the approx- 
imations above closes the resulting system of equations. The coefficient matrices, 
when written in the form of equation (7), are of order (N + 1) x (N + 1). The 
homogeneity of the boundary conditions, however, allows us to reduce the order of 
the lamda matrix by column operations to (N — 1) x (N — 1). The form of the 
boundary conditions are the same for both the collocation and the finite difference 
methods. They are given by 

v(±l) = 0 (18) 
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The appropriate column operations and substitutions can once again reduce the 
order of the coefficient matrices to (JV — 1) x (JV — 1). 

2.5 Eigenvalue and Eigenfunction Calculations 


Two methods are used to solve the eigenvalue problem, equation (7), in which 
the parameter is of order three. The LCM linearizes the problem and reduces 
it to a general eigenvalue problem. But the resulting companion matrix for the 
matrix polynomial is of higher order, 3 x ( N — 1), in the present case. The method, 
nevertheless, provides an approximation to the complete eigenvalue spectrum. The 
matrix factorization method features only a subset of the eigenvalue spectrum. The 
method, however, involves only matrices of order (JV - 1). The details of these two 
methods can be found in Bridges and Morris 10 . 

In general, the continuous part of the eigenvalue spectrum may be ignored 
when seeking only criterion for stability 18 . In order to numerically distinguish the 
discrete part from the continuous part of the spectrum, a transformation is used in 
conjunction with the matrix factorization method. The transformation is 


1 

(a f - a) 


(19) 


The lambda matrix then becomes 


D 3 (a) = C 0 6? + Cid 2 + C 2 d + C 3 (20) 

This transformation would insure that the eigenvalues of D 3 (o) in the vincinity of 
a f appear in the set of eigenvalues of the dominant solvent of D 3 (a). A solvent 
of D 3 (a) is said to be dominant if every^ eigenvalue in the set has an absolute 
value greater than all the eigenvalues of t) 3 (a) that are not in it. An algorithm 
developed by Dennis et al. 19 has been used to find the dominant solvent of the 
matrix polynomial. The desired eigenvalue can then be identified. The eigenvalues 
obtained from each of the previous methods may be further refined by the iterative 
method of Lancaster 20 . To compute single eigenvectors the inverse iteration method 
is used, 


D(a*) v* +1 = gv k (21) 

where g is a scaling factor and is taken as 1 in the calculation. The initial guess for 
the iteration is the complex unit vector. 

In this section, we have summarized three global approaches to solve the eigen- 
value problem generated by the spatial stability analysis of a free shear layer, in 
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which the eigenvalue appears in a nonlinear manner. The eigenvalues will be cal- 
culated using both the LCM and the matrix factorization methods. The results of 
our computations are presented in the next section. 

3. RESULTS AND DISCUSSION 

In this section, a comprehensive comparison of the different numerical methods 
is presented in terms of their accuracy and efficiency in solving the spatial, inviscid 
instability of a free mixing layer. 

The velocity profile for the basic free shear layer is chosen to be 

u (v) = |( 1 + tanli(y) ) (22) 

The physical domain is transformed to the domain [-1,1] on which Chebyshev poly- 
nomials are defined. Grosch and Orszag 21 studied spectral solutions of differential 
equations in both semi-infinite and infinite domains and used three types of trans- 
formations; These included a domain truncation, an algebraic mapping and an 
exponential mapping. They found that all of the three transformations are useful if 
the exact solution of the original differential equation decays exponentially fast as 
| y | — > oo, but fail if solutions oscillate out to infinity. Their results also showed 
that when the solution of a problem is smooth in the mapped domain algebraic map- 
pings are preferred over the other two. Since the solution of the Rayleigh equation 
in the regions of constant mean flow properties can be written as, 

as y — > ±oo (23) 


exp 


^fay 


a square-root transformation is used in the present analysis to avoid the singulari- 
ties at the end points of the transformed domain, 2 = ±1, which would arise if an 
exponential mapping was used 21 . As will be seen later, the square- root transforma- 
tion also enables the CC and the FD methods to predict accurately the continuous 
spectrum associated with the unboundedness even though the corresponding eigen- 
solutions are highly oscillatory at infinity. The transformation used is, 


2 = 


y 


where r is a scali 


( r 2 + f ) I/2 
factor. The metric dz/dy is 


m — 


(1 


5 ) 3/2 


(24) 


(25) 
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The scaling factor r controls the distribution of grid points. Increasing r decreases 
the number of grid points clustered around y = 0. Since the scaling factor determines 
the amount of the domain stretching, its optimum value, for which the solutions 
are most accurate, may depends on both the number of grid points used and the 
discretization scheme. However, the best grid distribution should be the same for 
a given problem, irrespective of the discretization scheme. Boyd 22 used a steepest 
descent method to predict the optimum choice of the mapping parameter in applying 
a Chebyshev polynomial approximation to a known, explicit function. Nevertheless, 
computing analytically the optimum mapping parameter in the application of global 
approximations to a differential equation is difficult. Some preliminary n um erical 
tests were performed using the matrix factorization method to evaluate the effects 
of domain stretchings on the prediction of the discrete eigenvalue spectrum . For 
= 0.2 and N — 17, the results are given in Table 1. The error, e, in each case is 
based upon the corresponding solutions from a shooting method and is determined 
by 


e 


, a 


a , 


a, 


(26) 


where a s is the eigenvalue calculated by the shooting method. This yields the value 


a;* = 0.38260 - iO. 22762. (27) 

In the shooting method the Rayleigh equation was integrated in the interval y € 
[—6,6] in 200 steps using a fourth-order, fixed step size Runge-Kutta procedure. 
Note that the solution decays exponentially at the far field and the center region 
around y = 0 is the region where there are large changes of the solution. To start 
a calculation, therefore, one can set r = 1. Since 


. m, (z = 0) = 1.0 (28) 

the region around z — 0, where there are large changes of flow properties, is not 
scaled for r — 1. Table 1 shows that the best scaling factors for the FD, the CC 
and the CT methods are 2.5, 2.0 and 2.5, respectively. Since the finite differencing 
is performed on equally spaced grids and the collocation points in the collocation 
method cluster at both ends of the computational domain, more domain stretching, 
or a bigger r , is needed for the finite differencing than for the collocation method. 
With N = 17, however, the dependence on the stretching parameters seems rather 
weak. The eigenvalues predicted by all the discretization methods are within 1% of 
the value obtained by using the shooting method. As will be shown later, this weak 
dependence quickly disappears as N increases. In the following calculations, the 
scaling factors used are the “optimum values” for each case unless otherwise noted. 

Tables 2 and 3 give the order of the coefficient matrices required by each method 
to obtain '10% and 1% of accuracy in the eigenvalue calculations, respectively. 
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Table 1 Predicted Eigenvalues. u> = 0.2, N — 17. 


r 

F D 

e 

C C 

e 

C T 

e 

0.5 

(0.38053,-0.14657) 

0.18 

(0.41107,-0.24395) 

0.073 

(0.36949,-0.18014) 

0.11 

1.0 

(0.35426,-0.22327) 

0.064 

(0.38119,-0.22371) 

0.009 

(0.37514,-0.23469) 

0.023 

1.5 

(0.37896,-0.23947) 

0.028 

(0.38171,-0.22870) 

0.003 

(0.38455,-0.22343) 

0.006 

2.0 

(0.38575,-0.22985) 

0.009 

(0.38290,-0.22834) 

0.002 

(0.38304,-0.22756) 

0.001 

2.5 

(0.38423,-0.22496) 

0.007 

(0.38313,-0.22931) 

0.004 

(0.38265,-0.22819) 

0.001 

3.0 

(0.38231,-0.22319) 

0.01 

(0.38559,-0.23188) 

0.012 

(0.38319,-0.22877) 

0.003 

4.0 

(0.38114,-0.22187) 

0.013 

(0.37169,-0.20992) 

0.047 

(0.38615,-0.23072) 

0.011 

5.0 

(0.38299,-0.22050) 

0.016 

(0.34586,-0.23982) 

0.087 

(0.38964,-0.23328) 

0.034 












It can be seen that the FD discretization predict the discete eigenvalue to the same 
order of accuracy as the CT and the CC methods do with the lowest order of the 
coefficient matrices. The two spectral methods performed almost equally well with 
proper choices of the mapping parameters. All the three discretization methods 
show rapid convergence at the low values of N. 

Table 2 Predicted Eigenvalues of less than 10% error, u = 0.2. 



N 

r 

O f 

— Oii 

e x 10 2 

FD 

5 

2.0 

0.38334 

0.19239 

7.9 

CC 

8 

1.5 

0.38870 

0.21940 

2.2 

CT 

7 

2.5 

0.41960 

0.20607 

9.6 

Predicted Eigenvalues of less than 1% error. 

u> = 0.2. 


N 

r 

Of y 

~0(i 

ex 10 2 

FD 

9 

2.5 

0.38364 

0.22812 

H ■ 

CC 

11 

1.5 

0.38264 

0.22364 

mSm 

CT 

10 

3.0 

0.38498 

0.22597 

0.6 


For the CC and the CT methods, the calculated eigenvalues converge from 10% to 
1% error by increasing the number of the approximation functions by about 30%. 
The reason is that in the mapped domain equation (23) becomes 


exp 


arz 


( l -* 2 ) 1 / 2 J 


as z 


±1 


(29) 


The solution is thus smooth in the mapped domain and the rapid convergence 
property of the global methods is retained. Comparisons of the rates of convergence 
for the various discretizations and solution techniques are given for the case to = 0.2 
in Table 4. All of the discretization methods show rapid (faster than algebraic) 
rates of convergence when the LCM is applied. Similarly, the MF method also gives 
the same rapid rate of convergence using the CT and the CC methods. As was 
expected, the finite difference discretization predicts the eigenvalues well but with 
lower rates of convergence. 

Note that the maximum computer time required in all the calculations in Table 
4 is less than a minute on a VAX 8550 machine. In practical applications, therefore, 
the effect of the mapping parameter can be effectively minimized by increasing the 
order of approximations, or N , without a significant increase in computer time. 
This is ^.lso evident in the following eigenfunction calculations. Figures 1 and 2 
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Figure 2. Eigenfunction for w = 0.2, JV = 23 using , Shooting Method; 

— — , CT; o , CC; x , FD. (a) Real Part, (b) Imaginary Part 







show the real and the imaginary parts of the calculated eigenfunctions for cases with 
u = 0.2 and iV=13 and 23, respectively. The error, e, of the corresponding calcu- 
lated eigenvalues was less than 1%. The agreement between the results is excellent 
even in the cases with 1V=13. For the discrete mode, all the methods become less 
sensitive to the choice of the scaling factors as N increases. While the spectral meth- 
ods performed better than the finite difference method in other calculations 6 ’ 23 , the 
finite difference formulation is competitive in the present application. Thus, the 
choice of solution methods appears to be problem- dependent. 

As can be seen from Table 5, the eigenvalues calculated by using the LCM 
argee to the fourth digits with those from the MF method for all the discretization 
methods. 


Table 5 Predicted Eigenvalues, lo = 0.2, N = 17. 



LCM 

MF 

r 

e x 10 2 

FD 

(0.38421,-0.22494) 

(0.38423,-0.22496) 

2.5 

0.7 

CC 

(0.38289,-0.22834) 

(0.38290,-0.22834)" 

2.0 

0.2 

CT 

(0.38304,-0.22755) 

(0.38304,-0.22756) 

2.0 

0.1 


The a / in the matrix factorization method are chosen such that 

(Xj = a 3 . (30) 

and the eigenvalues thus obtained are not refined by iterative methods. The a a , 
however, are not always known a priori for other flow conditions. For example, the 
realistic velocity profiles of mixing layers may be different from the one assumed 
here and, therefore, their eigenvalues may be different. It may thus be difficult to 
obtain converged eigensolutions using the shooting method. The eigenvalues and 
the corresponding e for other choices of a/ axe shown in Table 6. Both the FD 
and the CT methods give good results for up to 30% under-shoot of aif, for which 
the shooting method would have failed had the a / been used as an initial guess. 
Therefore, in conjunction with the FD and the CT methods, the MF method is far 
less sensitive to the choice of a f than the shooting method is to its initial guess. 
On the other hand, the Chebyshev collocation discretization is more sensitive to 
the a/. 

Figure 3 shows the growth rates of the spatially unstable modes of the free 
shear layer obtained by the various methods for N — 11. It was found that since 
the eigenfunctions decay more slowly in the far field as the frequency decreases, the 
mapping parameters selected for the mid frequency waves were not appropriate for 
the high and the low frequency wave calculations. The stretching parameters used 
were thus greater in the low frequency cases and smaller for the high frequencies. 
As was discussed earlier, the dependence quickly diminishes as N increases. 
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Table 6 . Eigenvalues and errors (e x lO 2 ) predicted using the matrix 
factorization method. 



a 

e x 10 2 

ccf - .(-0.3, 0.2), 1 f „ ' |X 10 2 - 20.0 

F D 

(0.376827,-0.224554) 

1.4 

c c 

(0.283691,-0.203911) 

22.0 


(0.377169,-0.220236) 

2.1 

a/ m (-0.228,0.132), 1 1 x 10 - 30 *° 

CX| 

F D 

(0.373056,-0.214566) 

3.6 

C C 

(0.232942,-0.263733) 

34.5 

C T 

(0.374547,-0.204189) 

in 

• 

in 
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One of the advantages of using the global approaches is that we obtain an 
approximation to the complete eigenvalue spectrum, while in the shooting method 
each eigenvalue has to be determined separately. In addition to a finite number 
of discrete values of a that satisfies the dispersion relation (3), there are also con- 
tinuous branches associated with the singularity of the Rayleigh equation. In the 
present stability calculations, for example, there are two continuous branches of 
the eigenvalue spectrum arising from the singularity, y c , of the Rayleigh equation, 
where 

a U(y c ) - to = 0. (31) 

For the hyperbolic-tangent velocity profile assumed here, these singular spectra are 

(i) ai = 0 ; a r > u (32) 

(ii) e R ; a r — > oo. (33) 

Figure 4a shows the calculated results with N — 17 and u>—0.2 using the FD and 
the LCM methods. The spectrum associated with the equations (32) can be clearly 
seen; however, the one associated with the equation (33) is cut out due to the 
magnitude of eigenvalues. 

Another continuous spectrum that can be observed in Figure 4 is associated 
with the bounded solutions of the asymptotic form of the Rayleigh equation in the 
far field, 

v" — a 2 v — 0. (34) 

since U" — > 0 as y — > droo. This is the continuous spectrum associated with the 
domain unboundedness and can be written as 

’ a r — 0, a; (E R (35) 

The corresponding eigenfunctions are purely perodic. As is shown in Figure 4b the 
numerical results predict better both of the continuous spectra, equation (32) and 
(35), with larger values of N. The finite approximations used here will converge to 
the solution of the equation as iV —» oo. Similar characteristics of the eigenvalue 
spectrum have been observed in plane Couette flow calculations 16 . The discrete 
part of the spectrum is associated with the convective instability. For u>=0.2, the 
plane Couette flow is stable and the discrete spectrum is empty. Case 24 obtained 
eigenfunctions corresponding to the eigenvalues in these singular branches by taking 
a Fourier transform with respect to x and a Laplace transform with respect to time 
of the linearized disturbance equations. The resulting equation is then solvable 
using a Green’s function method. The method is applicable to the current cases, 
but will not be given here. 
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Figure 4. Eigenvalue Spectra, w = 0.2. x , Eigenvalues Associated with the 
Unbounded Domain; a , Eigenvalues Associated with the Singularity 
of the Rayleigh Equation; • , Discrete Spectrum, (a) N = 17 (b) 
N = 76. ’ 







Figure 4 also shows that the presence of these continuous eigenvalue spectra 
may conceal or mask the discrete eigenvalues. The spurious solutions, however, 
are far away from the discrete spectrum in the complex wave speed plane. Figure 
5 shows this tendency clearly. The complex wave-speed of the discrete mode for 
oj—0.2 is (0.51845,-0.87404). Thus, the discrete spectrum can be better observed in 
the wave- speed plane. 

If a transformation that produces singularities at the boundaries of the trans- 
formed domain is used, the convergence property of the global approximations would 
no longer be retained. Figure 6 shows the approximated eigenvalue spectrum us- 
ing a hyperbolic-tangent transformation the CT method and the LCM method for 
lo = 0.2 and N = 27. The transformation is 

z = tanh (y). (36) 

Both the discrete and the continuous spectra are not well predicted even with the 
relative high order of approximation. 

Figures 7a, 7b and 7c show the eigenvalue spectrum with u> = 0.2 using the 
CC, the FD and the CT methods, respectively, and the square-root transformation. 
The spectra associated with the equation singularity, equations (32) and (33) and 
the discrete spectra are well predicted. As was discussed earlier, the convective 
instability described by the discrete spectrum is associated with the local vortic- 
ity distribution and is less sensitive to r as N increases. It can be observed from 
the present results that the same is true for the continuum due to the equation 
singularity. Despite the oscillatory nature of the corresponding eigenfunctions as 
y ~ * oo, the locations of the continuous spectra associated with the domain un- 
boundedness predicted by the CC and the FD also agree well with the analytic 
expression, equation (35). The square-root transformation used here seems to be a 
viable alternative to the three transformations tested by Grosch and Orszag 21 for 
either the FD or the CC methods. However, the continuum predicted by the CT 
method is very sensitive to the mapping parameter, r, even for the relatively high 
values of N. This may be due to the aliasing terms that are not included in the 
Chebyshev spectral tau methods. The eigenvalue spectrum is made complete with 
the inclusion of the continuous branches and an arbitrary initial disturbance can 
not be represented without knowing the complete eigenvalue spectrum. Therefore, 
a good approximation of the eigenvalue spectrum associated with the domain un- 
boundedness is important to the solution of the Rayleigh equation in an unbounded 
domain. 
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Figure 7. Eigenvalue Spectra for co — 0.2. * , N — 26, r = 2.0; a , JV = 
46, r = 2.0; o , TV = 76, r = 2.0; * , N = 76, r = 0.5. (a) Cheby- 
shev Collocation Method , (b) Finite Difference Method, (c) Chebyshev 
Tau Method. 










4. CONCLUSIONS 


Three discretization schemes, two spectral methods and a finite difference 
method have been applied to solve the spatial inviscid instability of a free mix- 
ing layer with a hyperbolic-tangent velocity profile. Calculated eigenfunctions for 
the discrete mode using these global approximations show good agreement with 
that using a conventional shooting procedure. For the same order of accuracy of 
the calculated eigenvalues, when compared to that using the shooting method, the 
finite difference discretization is more efficient than both the CKebyshev t'qu and tlie 
Chebyshev collocation methods. On the other hand, the Chebyshev tau method is 
more efficient than the Chebyshev collocation method. The finite difference method 
is also easier to formulate and code. All of the three discretization schemes result in 
rapid rates of convergence when the LCM is used. The matrix factorization is less 
sensitive to the atf than the shooting method is to its initial guess. The 0 ) / appears 
in the transformation that was used to identify the discrete eigenvalue. The LCM 
is preferred when the eigenvalue desired is not known a priori. The discrete part of 
the eigenvalue spectrum is very distinguishable in the complex wave speed plane. 

All of the discretization methods used here, the second-order finite difference, 
the Chebyshev tau and the Chebyshev collocation methods, are capable of pre- 
dicting accurately the discrete spectrum and the continuous spectrum associated 
with the singularity of the Rayleigh equation. The continuous spectrum associated 
with the unbounded domain can also be well predicted by the three method's, even 
though the Chebyshev tau predictions are somewhat more sensitive to the mapping 
parameter in the square-root transformation. The global eigensolution methods 
studied here may be applied very efficiently to obtain either an approximation to 
the complete eigenvalue spectrum or initial 1 gtiesses for a local shooting procedure 
for the discrete part of the spectrum. - ■ 
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